Supplementary information to report “How to measure response diversity”
Samuel R.P-J. Ross1, Owen L. Petchey, Takehiro Sasaki, David W. Armitage
27 January, 2023
rm(list = ls())
library(tidyverse)
library(patchwork)
library(ggtext)
library(mgcv)
library(gratia)
library(hypervolume)
library(here)
source(here("r/our_functions.r"))- Samuel R.P-J. Ross (Corresponding author), Integrative Community Ecology Unit, Okinawa Institute of Science & Technology Graduate University, Okinawa, Japan. E: s.ross.res[at]outlook.com
1 Introduction
This is the supplementary information accompanying the manuscript, How to measure response diversity by Samuel R.P-J. Ross, Owen L. Petchey, Takehiro Sasaki, and David W. Armitage.
The purpose of this document is to demonstrate how to measure the two dimensions of response diversity we propose in the manuscript, and their behavioural properties under a range of conditions (Section 2). We then apply these measures to a range of simulated examples (Section 3), and a reanalysis of a previous response diversity study using real data (Section 4). Finally, we provide description of data format for the results from our systematic review, and call those results (also provided as a .csv file).
2 Two dimensions of response diversity
A set of communities were created in which the species richness of the communities varied. Each of the \(i\) species in a community was give a trait value, in this case the first derivative of it’s response to an environmental variable (\(x_i\)).
## Create the communities
num_spp <- c(2, 4, 8, 16)
max_minus_min <- round(seq(1, 15, length = 45),3)
mean <- round(seq(-10, 10, length = 49),3)
rd_vals <- function(num_spp, max_minus_min, mean) {
seq(mean-max_minus_min/2, mean+max_minus_min/2, length = num_spp)
}
dd <- crossing(num_spp,
max_minus_min,
mean) %>%
rowwise() %>%
mutate(rd_vals = list(rd_vals(num_spp, max_minus_min, mean))) %>%
unnest(cols = rd_vals)The two dimensions (at least) of the response diversity of these values are:
To capture the range and spread of the data, we use a similarity-based diversity estimator described in detail in (Leinster & Cobbold 2012, Ecology, 93, 477-489). Mathematically, the similarity-based diversity of order \(q\), \(^qD^Z(p)\),is calculated as \(^qD^Z(p)=(\sum p_i(Zp)^{q-1}_i)^{\frac{1}{1-q}}\), where the exponent \(q\) defines the relative weighting of richness and evenness as is typical of Hill numbers. Here, \(q\) is fixed at zero, making diversity equivalent to species richness when \(Z=I\).The vector \(p\) contains the relative abundances for each species from \(1 \dots S\), and \(Z\) is an \(S \times S\) matrix encoding the similarity between the \(i^{th}\) and \(j^{th}\) species. Therefore, \((Zp)_i=\sum^S_{j=1}Z_{ij}p_j\). We can use the diversity, \(^qD^Z(p)\), as a metric capturing information on both richness and similarity in species’ responses (see Leinster & Cobbold 2012). There are many ways to calculate similarity, but here we use \(Z_{ij}=e^{-d_{ij}}\), where \(d_{ij}\) are the pairwise Euclidean distances (that is, dissimilarity) in performance-environment responses between two species. Thus, values of \(Z_{ij}\) can range from \(0\) as \(d_{ij} \to \infty\) to \(1\) when \(d_{ij}=0\). Further, cases where \(q \to 1\) and \(q \to \infty\) are defined as \(^1D^Z(p)=[(Zp)^{p_1}_1(Zp)^{p2}_1 \dots (Zp)^{p_S}_S]^{-1}\) and \(^{\infty}D^Z(p)=max(Zp)_i ^{-1}\), respectively.
To quantify the overlap with zero, we use \(\frac{range(x_i) - abs( abs(max(x_i)) - abs(min(x_i)) ) }{ range(x_i)}\). This quantity is zero when there is no overlap with zero, and is 1 when \(min(x_i) = -max(x_i)\).
We now demonstrate the properties of these indices under a range of conditions.
## calculate the diversity metrics
summary_stats <- dd %>%
group_by(num_spp, max_minus_min, mean) %>%
summarise(dissimilarity = resp_div(rd_vals,sign_sens = FALSE),
divergence = resp_div(rd_vals,sign_sens = TRUE)) %>%
pivot_longer(names_to = "variable", values_to = "value", 4:5)
## wrangle the data for easier plotting
all_data <- dd %>%
mutate(variable = "first_deriv_vals",
value = rd_vals) %>%
select(-rd_vals) %>%
bind_rows(summary_stats) %>%
mutate(variable = as_factor(variable))2.1 Communities varying only in spread
this_mean <- 0
this_num_spp <- 8
all_data %>%
filter(mean == this_mean,
num_spp == this_num_spp) %>%
ggplot(aes(x = max_minus_min, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Range of first derivative values in community") 
Figure 2.1: Communities varying only in the range of trait values. Trait values overlap zero. but the extent of overlap is constant, number of species is constant, and mean trait value is constant. Each column of dots corresponds to one community. Top panel shows the values of the first derivative of each of the species in the community. Middle panel is the measure of dissimilarity. Bottom panel is the measure of divergence.
this_mean <- 2.5
this_num_spp <- 8
all_data %>%
filter(mean == this_mean,
num_spp == this_num_spp,
max_minus_min < 5) %>%
ggplot(aes(x = max_minus_min, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Range of first derivative values in community")
Figure 2.2: Communities varying only in the range of trait values, with no overlap of trait values with zero. Extent of overlap of trait values with zero is constant, number of species is constant and mean trait value is constant. All else is as in the previous figure.
2.2 Communities varying only in overlap
Here we look at communities that only vary in the overlap of trait values with zero, and that do not differ in range. The measure of the overlap with zero is zero while there is not overlap, then increases to a maximum when the mean is zero, and then decreases again (Fig. 2.3.
this_max_minus_min <- 5.455
this_num_spp <- 8
all_data %>%
filter(max_minus_min == this_max_minus_min,
num_spp == this_num_spp) %>%
ggplot(aes(x = mean, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Mean of first derivative values in community") 
Figure 2.3: Communities varying only in the mean trait value, but not in range of trait values or the number of species.
2.3 Communities varying in both spread and overlap
this_mean <- 2.5
this_num_spp <- 8
all_data %>%
filter(mean == this_mean,
num_spp == this_num_spp) %>%
ggplot(aes(x = max_minus_min, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Range of first derivative values in community")
Figure 2.4: Communities varying both the range and the amount of overlap, but with constant mean trait value.
2.4 Varying number of species
When only number of species varies, dissimilarity increases.
this_mean <- 0
this_num_spp <- 8
this_max_minus_min <- 5.455
all_data %>%
filter(mean == this_mean,
max_minus_min == this_max_minus_min) %>%
ggplot(aes(x = num_spp, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Number of species") 
Figure 2.5: Communities with different numbers of species, but same range of trait values, same amount of overlap with zero, and same mean trait values.
2.5 Redundant species additions
When species are added with identical responses to species already in the community, dissimilarity does not change.That is, redundancy does not inflate response diversity (another desirable property)
num_spp <- c(2, 4, 8, 16,3)
redundant<-data.frame('S'=num_spp %>% rep(times=c(2,4,8,16,3)),'first_deriv_vals'=NA,'dissimilarity'=NA,'divergence'=NA)
redundant$first_deriv_vals[c(1:32)]<-c(1,-1)
redundant$first_deriv_vals[33]<-1
for(i in 1:length(redundant$S)){
first_deriv_vals<-c(1,-1) %>% rep(redundant$S[i]*0.5)
redundant$dissimilarity[i]<-first_deriv_vals %>% resp_div(sign_sens = F)
redundant$divergence[i]<-first_deriv_vals %>% resp_div(sign_sens = T)
}
redundant<-redundant %>%
pivot_longer(names_to = "variable", values_to = "value", 2:4)
redundant$variable<-redundant$variable %>%
as.character() %>% parse_factor(levels = c('first_deriv_vals','dissimilarity','divergence'))
redundant %>%
ggplot(aes(x = S, y = value)) +
geom_point() +
geom_hline(yintercept = 0, lty="dashed") +
facet_wrap( ~ variable, scales = "free_y", ncol = 1) +
xlab("Number of species") 
Figure 2.6: Communities with different numbers of species, but identical trait values.
3 Producing Figures 2-3 based on simulated data
This section will produce Figures 2-3 from the manuscript. The figures are inspired by a range of examples from the literature (see main text), but are based on simulated data.
First, we’ll generate a series of simulated communities that have high, medium, or low response diversity to an environmental variable [x], based on correlated or uncorrelated responses, and responding either linearly or nonlinearly.
# generate an environmental variable (x), as a simple sequence of numbers
x = seq(0, 1, length.out = 100)
# generate species with different responses to the environment
# Linear responses:
y1a <- 0.8-1*x; y2a <- 0.75-1* x; y3a <-0.85-1*x # responses are similar in direction and magnitude (low response diversity)
y1b <- 0.5-0.2*x; y2b <- 0.75-1*x; y3b <-1-2.5*x # responses are similar in direction but not magnitude (medium response diversity)
y1c <- 0.6-0.75*x; y2c <- rev(y1c); y3c <- rep(0.35, length(x)) # responses are not similar in direction or magnitude (high response diversity)
# Nonlinear responses:
y1d <- range01(dnorm(x,0.5,0.3)); y2d <- range01(dnorm(x,0.475,0.3)); y3d <- range01(dnorm(x,0.525,0.3)) # responses are similar in direction and magnitude (low response diversity)
y1e <- range01(dnorm(x,0.5,0.3)); y2e <- range01(dnorm(x,0.2,0.3)); y3e <- range01(dnorm(x,0.8,0.3)) # responses are similar in direction but not magnitude (medium response diversity)
y1f <- exp(-3*x); y2f <- rev(y1f); y3f <- 1-(y1f+y2f) # responses are not similar in direction or magnitude (high response diversity)We’ll model these species-environment relationships using generalized additive models (GAMs).
# model individual species~environment relationships using GAM
m1a <- gam(y1a~s(x)); m2a <- gam(y2a~s(x)); m3a <- gam(y3a~s(x))
m1b <- gam(y1b~s(x)); m2b <- gam(y2b~s(x)); m3b <- gam(y3b~s(x))
m1c <- gam(y1c~s(x)); m2c <- gam(y2c~s(x)); m3c <- gam(y3c~s(x))
m1d <- gam(y1d~s(x)); m2d <- gam(y2d~s(x)); m3d <- gam(y3d~s(x))
m1e <- gam(y1e~s(x)); m2e <- gam(y2e~s(x)); m3e <- gam(y3e~s(x))
m1f <- gam(y1f~s(x)); m2f <- gam(y2f~s(x)); m3f <- gam(y3f~s(x))Then we calculate the first derivatives of these relationships.
# calculate first derivatives from models
d1a <- derivatives(m1a); d2a <- derivatives(m2a); d3a <- derivatives(m3a)
d1b <- derivatives(m1b); d2b <- derivatives(m2b); d3b <- derivatives(m3b)
d1c <- derivatives(m1c); d2c <- derivatives(m2c); d3c <- derivatives(m3c)
d1d <- derivatives(m1d); d2d <- derivatives(m2d); d3d <- derivatives(m3d)
d1e <- derivatives(m1e); d2e <- derivatives(m2e); d3e <- derivatives(m3e)
d1f <- derivatives(m1f); d2f <- derivatives(m2f); d3f <- derivatives(m3f)Then we’ll measure response diversity using the approach above.
# calculate response diversity for model a
rdiv_a<-data.frame('data' = d1a$data, 'sp.1' = d1a$derivative, 'sp.2' = d2a$derivative, 'sp.3' = d3a$derivative)
rdiv_a$rdiv<-apply(rdiv_a[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_a$sign<-apply(rdiv_a[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_a$Med<-median(rdiv_a$rdiv)
# calculate response diversity for model b
rdiv_b<-data.frame('data' = d1b$data, 'sp.1' = d1b$derivative, 'sp.2' = d2b$derivative, 'sp.3' = d3b$derivative)
rdiv_b$rdiv<-apply(rdiv_b[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_b$sign<-apply(rdiv_b[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_b$Med<-median(rdiv_b$rdiv)
# calculate response diversity for model c
rdiv_c<-data.frame('data' = d1c$data, 'sp.1' = d1c$derivative, 'sp.2' = d2c$derivative, 'sp.3' = d3c$derivative)
rdiv_c$rdiv<-apply(rdiv_c[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_c$sign<-apply(rdiv_c[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_c$Med<-median(rdiv_c$rdiv)
# calculate response diversity for model d
rdiv_d<-data.frame('data' = d1d$data, 'sp.1' = d1d$derivative, 'sp.2' = d2d$derivative, 'sp.3' = d3d$derivative)
rdiv_d$rdiv<-apply(rdiv_d[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_d$sign<-apply(rdiv_d[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_d$Med<-median(rdiv_d$rdiv)
# calculate response diversity for model e
rdiv_e<-data.frame('data' = d1e$data, 'sp.1' = d1e$derivative, 'sp.2' = d2e$derivative, 'sp.3' = d3e$derivative)
rdiv_e$rdiv<-apply(rdiv_e[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_e$sign<-apply(rdiv_e[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_e$Med<-median(rdiv_e$rdiv)
# calculate response diversity for model f
rdiv_f<-data.frame('data' = d1f$data, 'sp.1' = d1f$derivative, 'sp.2' = d2f$derivative, 'sp.3' = d3f$derivative)
rdiv_f$rdiv<-apply(rdiv_f[,c(2:4)], 1, resp_div, sign_sens = F)
rdiv_f$sign<-apply(rdiv_f[,c(2:4)], 1, resp_div, sign_sens = T)
rdiv_f$Med<-median(rdiv_f$rdiv)Now we have all the data required to produce the figures from the manuscript. First, we’ll produce Figure 2, which looks at linear performance-environment relationships:
# Y value range
ymax<-max(y1a,y2a,y3a,y1b,y2b,y3b,y1c,y2c,y3c);ymin<-min(y1a,y2a,y3a,y1b,y2b,y3b,y1c,y2c,y3c)
dmax<-max(d1a$derivative,d2a$derivative,d3a$derivative,d1b$derivative,d2b$derivative,d3b$derivative,d1c$derivative,d2c$derivative,d3c$derivative)
dmin<-min(d1a$derivative,d2a$derivative,d3a$derivative,d1b$derivative,d2b$derivative,d3b$derivative,d1c$derivative,d2c$derivative,d3c$derivative)
rmax<-max(rdiv_a$rdiv,rdiv_b$rdiv,rdiv_c$rdiv);rmin<-1
smax<-1;smin<-0
# panel a
Fig_2a<-ggplot(data = NULL,mapping = aes(x = x,y = y1a)) +
theme_classic(base_size = 14) +
labs(x = "Spring frost",y = "Chlorophyll production",title = "Low response diversity",tag = "a)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2a),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3a),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel b
Fig_2b<-ggplot(data = NULL,mapping = aes(x = x,y = y1b)) +
theme_classic(base_size = 14) +
labs(x = "Urbanisation",y = "Macroinvertebrate abundance",title = "Medium response diversity",tag = "b)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2b),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3b),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel c
Fig_2c<-ggplot(data = NULL,mapping = aes(x = x,y = y1c)) +
theme_classic(base_size = 14) +
labs(x = "Winter temperature",y = "Bee mass loss",title = "High response diversity",tag = "c)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2c),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3c),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel d
Fig_2d<-ggplot(data = NULL,mapping = aes(x = d1a$data,y = d1a$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Spring frost",y = "Derivative (chlorophyll product.)",tag = "d)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2a$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3a$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel e
Fig_2e<-ggplot(data = NULL,mapping = aes(x = d1b$data,y = d1b$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Urbanisation",y = "Derivative (macroinvert. abundance)",tag = "e)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2b$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3b$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel f
Fig_2f<-ggplot(data = NULL,mapping = aes(x = d1c$data,y = d1c$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Winter temperature",y = "Derivative (bee mass loss)",tag = "f)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2c$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3c$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel g
Fig_2ga<-ggplot(data = rdiv_a,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "g)") +
geom_hline(yintercept = rdiv_a$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = rdiv_a$Med+0.25),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_a$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_2gb<-ggplot(data = rdiv_a,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Spring frost",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_2gb<-Fig_2gb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# panel h
Fig_2ha<-ggplot(data = rdiv_b,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "h)") +
geom_hline(yintercept = rdiv_b$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = rdiv_b$Med-0.25),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_b$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_2hb<-ggplot(data = rdiv_b,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Urbanisation",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_2hb<-Fig_2hb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# panel i
Fig_2ia<-ggplot(data = rdiv_c,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "i)") +
geom_hline(yintercept = rdiv_c$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = rdiv_c$Med+0.25),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_c$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_2ib<-ggplot(data = rdiv_c,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Winter temperature",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_2ib<-Fig_2ib + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
(Fig_2a + Fig_2b + Fig_2c)/(Fig_2d + Fig_2e + Fig_2f)/(Fig_2ga + Fig_2ha + Fig_2ia)/(Fig_2gb + Fig_2hb + Fig_2ib) + plot_layout(heights = c(3,3,3,1))
Next, we’ll produce Figure 3 to look at nonlinear performance-environment relationships:
# Y value range
ymax<-max(y1d,y2d,y3d,y1e,y2e,y3e,y1f,y2f,y3f);ymin<-min(y1d,y2d,y3d,y1e,y2e,y3e,y1f,y2f,y3f)
dmax<-max(d1d$derivative,d2d$derivative,d3d$derivative,d1e$derivative,d2e$derivative,d3e$derivative,d1f$derivative,d2f$derivative,d3f$derivative)
dmin<-min(d1d$derivative,d2d$derivative,d3d$derivative,d1e$derivative,d2e$derivative,d3e$derivative,d1f$derivative,d2f$derivative,d3f$derivative)
rmax<-max(rdiv_d$rdiv,rdiv_e$rdiv,rdiv_f$rdiv);rmin<-1
smax<-1;smin<-0
# panel a
Fig_3a<-ggplot(data = NULL,mapping = aes(x = x,y = y1d)) +
theme_classic(base_size = 14) +
labs(x = "Flowering plant abundance",y = "Flowering plant rate of increase",title = "Low response diversity",tag = "a)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2d),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3d),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel b
Fig_3b<-ggplot(data = NULL,mapping = aes(x = x,y = y1e)) +
theme_classic(base_size = 14) +
labs(x = "Grazing pressure",y = "Seedling establishment",title = "Medium response diversity",tag = "b)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2e),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3e),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel c
Fig_3c<-ggplot(data = NULL,mapping = aes(x = x,y = y1f)) +
theme_classic(base_size = 14) +
labs(x = "Fire frequency",y = "Tree density",title = "High response diversity",tag = "c)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = y2f),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = y3f),
col = Isfahan1[[1]][6]) +
lims(y = c(ymin,ymax))
# panel d
Fig_3d<-ggplot(data = NULL,mapping = aes(x = d1d$data,y = d1d$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Flowering plant abundance",y = "Derivative (rate of increase)",tag = "d)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2d$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3d$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel e
Fig_3e<-ggplot(data = NULL,mapping = aes(x = d1e$data,y = d1e$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Grazing pressure",y = "Derivative (seedling establish.)",tag = "e)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2e$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3e$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel f
Fig_3f<-ggplot(data = NULL,mapping = aes(x = d1f$data,y = d1f$derivative)) +
theme_classic(base_size = 14) +
labs(x = "Fire frequency",y = "Derivative (tree density)",tag = "f)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][2]) +
geom_line(aes(y = d2f$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3f$derivative),
col = Isfahan1[[1]][6]) +
lims(y = c(dmin,dmax))
# panel g
Fig_3ga<-ggplot(data = rdiv_d,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "g)") +
geom_hline(yintercept = rdiv_d$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = rdiv_d$Med+1.25),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_d$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_3gb<-ggplot(data = rdiv_d,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Flowering plant abundance",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_3gb<-Fig_3gb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# panel h
Fig_3ha<-ggplot(data = rdiv_e,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "h)") +
geom_hline(yintercept = rdiv_e$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = rdiv_e$Med-0.75),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_e$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_3hb<-ggplot(data = rdiv_e,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Grazing pressure",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_3hb<-Fig_3hb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# panel i
Fig_3ia<-ggplot(data = rdiv_f,mapping = aes(x = data,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "i)") +
geom_hline(yintercept = rdiv_f$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 0.85,
y = 1.25),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_f$Med,digits = 2)))) +
lims(y = c(rmin,rmax))
Fig_3ib<-ggplot(data = rdiv_f,mapping = aes(x = data,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Fire frequency",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_3ib<-Fig_3ib + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
(Fig_3a + Fig_3b + Fig_3c)/(Fig_3d + Fig_3e + Fig_3f)/(Fig_3ga + Fig_3ha + Fig_3ia)/(Fig_3gb + Fig_3hb + Fig_3ib) + plot_layout(heights = c(3,3,3,1))
Demonstrate the relationship between the dissimilarity and divergence metrics:
# linear
cor(rdiv_a$rdiv,rdiv_a$sign); plot(rdiv_a$rdiv,rdiv_a$sign)## [1] NA
cor(rdiv_b$rdiv,rdiv_b$sign); plot(rdiv_b$rdiv,rdiv_b$sign)## [1] NA
cor(rdiv_c$rdiv,rdiv_c$sign); plot(rdiv_c$rdiv,rdiv_c$sign)## [1] 0.01887203
# nonlinear
cor(rdiv_d$rdiv,rdiv_d$sign); plot(rdiv_d$rdiv,rdiv_d$sign)## [1] 0.3574677
cor(rdiv_e$rdiv,rdiv_e$sign); plot(rdiv_e$rdiv,rdiv_e$sign)## [1] 0.7047028
cor(rdiv_f$rdiv,rdiv_f$sign); plot(rdiv_f$rdiv,rdiv_f$sign)## [1] 0.745555
Note that in the linear cases, response diversity metrics are only weakly correlated, whereas for nonlinear examples, the relationship between the two metrics is stronger (particularly as average divergence increases).
4 Leary & Petchey (2009) reanalysis
This section reanalyses the data from Leary & Petchey (2009, Journal of Animal Ecology) using the framework we present above.
4.1 Calculating response diversity
rm(list = ls())
library(tidyverse)
library(patchwork)
library(ggtext)
library(mgcv)
library(gratia)
library(hypervolume)
library(here)
source(here("r/our_functions.r"))#rm(list = ls())
r_vals <- read_csv(here("data/Leary_and_Petchey_r_values.csv")) %>%
mutate(parameter = "r")
K_vals <- read_csv(here("data/Leary_and_Petchey_K_values.csv")) %>%
mutate(x = log10(x), parameter = "log10_K")
rK_vals <- rbind(r_vals, K_vals) %>%
rename(value = x)
#rm(r_vals, K_vals)The reanalysis is based on the range of the first derivatives and including information on whether those first derivatives span zero (and by how much).
# get subsets of the data for measuring response diversity
sp1_r<-rK_vals[rK_vals$species %in% 'tetrahymena' & rK_vals$parameter %in% 'r',]
sp2_r<-rK_vals[rK_vals$species %in% 'loxocephallus' & rK_vals$parameter %in% 'r',]
sp3_r<-rK_vals[rK_vals$species %in% 'colpidium' & rK_vals$parameter %in% 'r',]
sp4_r<-rK_vals[rK_vals$species %in% 'paramecium' & rK_vals$parameter %in% 'r',]
sp1_K<-rK_vals[rK_vals$species %in% 'tetrahymena' & rK_vals$parameter %in% 'log10_K',]
sp2_K<-rK_vals[rK_vals$species %in% 'loxocephallus' & rK_vals$parameter %in% 'log10_K',]
sp3_K<-rK_vals[rK_vals$species %in% 'colpidium' & rK_vals$parameter %in% 'log10_K',]
sp4_K<-rK_vals[rK_vals$species %in% 'paramecium' & rK_vals$parameter %in% 'log10_K',]
# using k=5 to avoid the error when the predictor variable has fewer levels than k in model spec. <https://stackoverflow.com/a/62914481>
m1_r<-gam(data = sp1_r,value~s(temperature,k=5))
m2_r<-gam(data = sp2_r,value~s(temperature,k=5))
m3_r<-gam(data = sp3_r,value~s(temperature,k=5))
m4_r<-gam(data = sp4_r,value~s(temperature,k=5))
m1_K<-gam(data = sp1_K,value~s(temperature,k=5))
m2_K<-gam(data = sp2_K,value~s(temperature,k=5))
m3_K<-gam(data = sp3_K,value~s(temperature,k=5))
m4_K<-gam(data = sp4_K,value~s(temperature,k=5))
## predict
new_data <- data.frame(temperature = seq(18, 26, by= 0.02))
#new_data <- data.frame(temperature = seq(18, 26, by= 1))
p1_K <- predict(m1_K, newdata = new_data)
# get first derivatives
d1_r<-derivatives(m1_r, newdata = new_data)
d2_r<-derivatives(m2_r, newdata = new_data)
d3_r<-derivatives(m3_r, newdata = new_data)
d4_r<-derivatives(m4_r, newdata = new_data)
d1_K<-derivatives(m1_K, newdata = new_data)
d2_K<-derivatives(m2_K, newdata = new_data)
d3_K<-derivatives(m3_K, newdata = new_data)
d4_K<-derivatives(m4_K, newdata = new_data)
# calculate response diversity for intrinsic rate of increase (r)
rdiv_r<-data.frame('temperature' = d1_r$data,
'tetrahymena' = d1_r$derivative,
'loxocephallus' = d2_r$derivative,
'colpidium' = d3_r$derivative, 'paramecium' = d4_r$derivative)
rdiv_r$rdiv<-apply(rdiv_r[,c(2:5)], 1, resp_div, sign_sens = F)
rdiv_r$sign<-apply(rdiv_r[,c(2:5)], 1, resp_div, sign_sens = T)
rdiv_r$Med<-median(rdiv_r$rdiv)
# calculate response diversity for carrying capacity (K)
rdiv_K<-data.frame('temperature' = d1_K$data,
'tetrahymena' = d1_K$derivative,
'loxocephallus' = d2_K$derivative,
'colpidium' = d3_K$derivative,
'paramecium' = d4_K$derivative)
rdiv_K$rdiv<-apply(rdiv_K[,c(2:5)], 1, resp_div, sign_sens = F)
rdiv_K$sign<-apply(rdiv_K[,c(2:5)], 1, resp_div, sign_sens = T)
rdiv_K$Med<-median(rdiv_K$rdiv)To plot the output of the gams, we’ll use the predict function. This needs a new tibble with the output of predict:
# generate new tibble with data from predict() for plotting gams ## r
gam_data_r <- tibble(temperature = seq(from=min(sp1_r$temperature),to=max(sp1_r$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m1_r, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(gam_data_r)[c(2:3)]<-c("Sp1","Sp1_se") # rename to species 1
temp<-tibble(temperature = seq(from=min(sp1_r$temperature),to=max(sp1_r$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m2_r, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp2","Sp2_se") # rename to species 2
gam_data_r<-gam_data_r %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
temp<-tibble(temperature = seq(from=min(sp1_r$temperature),to=max(sp1_r$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m3_r, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp3","Sp3_se") # rename to species 2
gam_data_r<-gam_data_r %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
temp<-tibble(temperature = seq(from=min(sp1_r$temperature),to=max(sp1_r$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m4_r, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp4","Sp4_se") # rename to species 2
gam_data_r<-gam_data_r %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
# add cols for raw data from Leary & Petchey
gam_data_r$sp1_raw<-NA
gam_data_r$sp2_raw<-NA
gam_data_r$sp3_raw<-NA
gam_data_r$sp4_raw<-NA
gam_data_r<-gam_data_r %>% bind_rows(as_tibble(
sp1_r %>% group_by(temperature) %>% summarise(sp1_raw = value)
))
sp2_r <- sp2_r %>% group_by(temperature) %>% summarise(value = value)
sp3_r <- sp3_r %>% group_by(temperature) %>% summarise(value = value)
sp4_r <- sp4_r %>% group_by(temperature) %>% summarise(value = value)
for (i in 1:nrow(sp2_r)) {
gam_data_r$sp2_raw[1000+i]<-sp2_r$value[i]
gam_data_r$sp3_raw[1000+i]<-sp3_r$value[i]
gam_data_r$sp4_raw[1000+i]<-sp4_r$value[i]
}
# generate new tibble with data from predict() for plotting gams ## K
gam_data_K <- tibble(temperature = seq(from=min(sp1_K$temperature),to=max(sp1_K$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m1_K, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(gam_data_K)[c(2:3)]<-c("Sp1","Sp1_se") # rename to species 1
temp<-tibble(temperature = seq(from=min(sp1_K$temperature),to=max(sp1_K$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m2_K, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp2","Sp2_se") # rename to species 2
gam_data_K<-gam_data_K %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
temp<-tibble(temperature = seq(from=min(sp1_K$temperature),to=max(sp1_K$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m3_K, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp3","Sp3_se") # rename to species 2
gam_data_K<-gam_data_K %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
temp<-tibble(temperature = seq(from=min(sp1_K$temperature),to=max(sp1_K$temperature),length.out = 1000)
) %>%
bind_cols(as_tibble(
predict(m4_K, newdata = ., se.fit=TRUE)
)) # get gam data
colnames(temp)[c(2:3)]<-c("Sp4","Sp4_se") # rename to species 2
gam_data_K<-gam_data_K %>% bind_cols(temp[c(2,3)]); rm(temp) # bind on sp 2 and remove temporary tibble
# add cols for raw data from Leary & Petchey
gam_data_K$sp1_raw<-NA
gam_data_K$sp2_raw<-NA
gam_data_K$sp3_raw<-NA
gam_data_K$sp4_raw<-NA
gam_data_K<-gam_data_K %>% bind_rows(as_tibble(
sp1_K %>% group_by(temperature) %>% summarise(sp1_raw = value)
))
sp2_K <- sp2_K %>% group_by(temperature) %>% summarise(value = value)
sp3_K <- sp3_K %>% group_by(temperature) %>% summarise(value = value)
sp4_K <- sp4_K %>% group_by(temperature) %>% summarise(value = value)
for (i in 1:nrow(sp2_K)) {
gam_data_K$sp2_raw[1000+i]<-sp2_K$value[i]
gam_data_K$sp3_raw[1000+i]<-sp3_K$value[i]
gam_data_K$sp4_raw[1000+i]<-sp4_K$value[i]
}We can then plot a figure similar to Figures 2 and 3, but instead based on the real data from Leary & Petchey (2009). This is Figure 4 (or Box Figure 1) from the manuscript.
smax<-1;smin<-0
# panel a
Fig_4a<-gam_data_r %>%
ggplot(aes(x=temperature, y=Sp1)) +
theme_classic(base_size = 14) +
labs(x = "Temperature",y = "Intrinsic rate of increase (r)",tag = "a)") +
geom_point(aes(y = sp1_raw),col = Isfahan1[[1]][5],size = 2,shape = 1) +
geom_point(aes(y = sp2_raw),col = Isfahan1[[1]][4],size = 2,shape = 1) +
geom_point(aes(y = sp3_raw),col = Isfahan1[[1]][3],size = 2,shape = 1) +
geom_point(aes(y = sp4_raw),col = Isfahan1[[1]][2],size = 2,shape = 1) +
geom_ribbon(aes(ymin=Sp1-Sp1_se, ymax=Sp1+Sp1_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp2-Sp2_se, ymax=Sp2+Sp2_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp3-Sp3_se, ymax=Sp3+Sp3_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp4-Sp4_se, ymax=Sp4+Sp4_se),alpha=.1) +
geom_line(col = Isfahan1[[1]][5]) +
geom_line(aes(y = Sp2), col = Isfahan1[[1]][4]) +
geom_line(aes(y = Sp3), col = Isfahan1[[1]][3]) +
geom_line(aes(y = Sp4), col = Isfahan1[[1]][2])
# panel b
Fig_4b<-gam_data_K %>%
ggplot(aes(x=temperature, y=Sp1)) +
theme_classic(base_size = 14) +
labs(x = "Temperature",y = "log10 carrying capacity (K)",tag = "b)") +
geom_point(aes(y = sp1_raw),col = Isfahan1[[1]][5],size = 2,shape = 1) +
geom_point(aes(y = sp2_raw),col = Isfahan1[[1]][4],size = 2,shape = 1) +
geom_point(aes(y = sp3_raw),col = Isfahan1[[1]][3],size = 2,shape = 1) +
geom_point(aes(y = sp4_raw),col = Isfahan1[[1]][2],size = 2,shape = 1) +
geom_ribbon(aes(ymin=Sp1-Sp1_se, ymax=Sp1+Sp1_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp2-Sp2_se, ymax=Sp2+Sp2_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp3-Sp3_se, ymax=Sp3+Sp3_se),alpha=.1) +
geom_ribbon(aes(ymin=Sp4-Sp4_se, ymax=Sp4+Sp4_se),alpha=.1) +
geom_line(col = Isfahan1[[1]][5]) +
geom_line(aes(y = Sp2), col = Isfahan1[[1]][4]) +
geom_line(aes(y = Sp3), col = Isfahan1[[1]][3]) +
geom_line(aes(y = Sp4), col = Isfahan1[[1]][2])
# panel c
Fig_4c<-ggplot(data = d1_r,mapping = aes(x = data,y = derivative)) +
theme_classic(base_size = 14) +
labs(x = "Temperature",y = "Derivative (rate of increase)",tag = "c)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][5]) +
geom_line(aes(y = d2_r$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3_r$derivative),
col = Isfahan1[[1]][3]) +
geom_line(aes(y = d4_r$derivative),
col = Isfahan1[[1]][2])
# panel d
Fig_4d<-ggplot(data = d1_K,mapping = aes(x = data,y = derivative)) +
theme_classic(base_size = 14) +
labs(x = "Temperature",y = "Derivative (carrying capacity)",tag = "d)") +
geom_hline(yintercept = 0,
lty=2) +
geom_line(col = Isfahan1[[1]][5]) +
geom_line(aes(y = d2_K$derivative),
col = Isfahan1[[1]][4]) +
geom_line(aes(y = d3_K$derivative),
col = Isfahan1[[1]][3]) +
geom_line(aes(y = d4_K$derivative),
col = Isfahan1[[1]][2])
# panel e
Fig_4ea<-ggplot(data = rdiv_r,mapping = aes(x = temperature,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "e)") +
geom_hline(yintercept = rdiv_r$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 19,
y = rdiv_r$Med+0.05),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_r$Med,digits = 2))))
Fig_4eb<-ggplot(data = rdiv_r,mapping = aes(x = temperature,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Temperature",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_4eb<-Fig_4eb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# panel f
Fig_4fa<-ggplot(data = rdiv_K,mapping = aes(x = temperature,y = rdiv)) +
theme_classic(base_size = 14) +
labs(x = NULL,y = "dissimilarity (derivatives)",tag = "f)") +
geom_hline(yintercept = rdiv_K$Med,lty=2) +
geom_line() +
geom_richtext(mapping = aes(x = 19,
y = rdiv_K$Med+0.05),
size=4.5,
label.color = NA,
label = paste("RDiv^Med =",paste0(round(rdiv_K$Med,digits = 2))))
Fig_4fb<-ggplot(data = rdiv_K,mapping = aes(x = temperature,y = sign)) +
theme_bw(base_size = 14)+
labs(x = "Temperature",y = "Divergence") +
geom_line() +
lims(y = c(smin,smax))
Fig_4fb<-Fig_4fb + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
(Fig_4a + Fig_4b)/(Fig_4c + Fig_4d)/(Fig_4ea + Fig_4fa)/(Fig_4eb + Fig_4fb) + plot_layout(heights = c(3,3,3,1))
Quick check to see if r and K are correlated:
cor(r_vals$x,K_vals$x,method = "spearman")## [1] 0.1299806Low Spearman coefficient (0.13).
4.2 Comparing response diversity measures
Here we read in the data that was created in the benchmarking response diversity project.
## read in the community level data
comm_data <- read_csv(here("data/all_comm_data.csv")) ## Have a look at relationships between correlation based measures and derivative based measures:
ggplot(comm_data) +
geom_point(aes(x = cor_respdiv, y=deriv_nss_respdiv,
col = environment_values, shape = composition)) +
facet_wrap(vars(trait_name))
ggplot(comm_data) +
geom_point(aes(x = cor_respdiv, y=deriv_ss_respdiv,
col = environment_values, shape = composition)) +
facet_wrap(vars(trait_name))We’ve looked at cases where the correlation based response diversity is either very positive or very negative, and the sign sensitive derivative is respectively either very high or very low (the opposite of what one would expect). We think this occurs mainly due to the correlation being calculated only over the temperature range experienced by the community, whereas the gam from which the derivatives are derived is fit to the entire temperature range, and then only the part of the fit for the temperature range experienced is used (even though that can be affected by the data points outside that range in the gam fit). So even before we calculate response diversity, we are making choices that affect the outcome (i.e. value of response diversity).
## check some of the positive correlations but sign difference
## K, 24, LT. Correlation is about 0.5, divergence is about 0.18
## plot the rs
LT_24_K <- K_vals %>%
group_by(species, temperature) %>%
summarise(x = mean(x)) %>%
filter(temperature > 20 , species %in% c("loxocephallus", "tetrahymena")) %>%
#select(-jar) %>%
pivot_wider(names_from = species, values_from = x)
ggplot(LT_24_K) +
geom_point(aes(x = tetrahymena, y = loxocephallus, col = temperature))
cor(LT_24_K$tetrahymena, LT_24_K$loxocephallus) ## not exactly the same as in Leary and Petchey (i.e. in comm_stats), but this is because Leary and Petchey used a specific method for calculating the correlations, which is not the same as this method.
## E.g. Petchey & Leary calculated a correlation for every possible pairing of the jars, and took the mean of that. And for carrying capacity, correlations were calculated on the raw biomass, and not the log10 biomass
LT_24_r <- r_vals %>%
filter(temperature > 20, species %in% c("loxocephallus", "tetrahymena")) %>%
select(-jar) %>%
pivot_wider(names_from = species, values_from = x)
ggplot(LT_24_r) +
geom_point(aes(x = tetrahymena, y = loxocephallus))
cor(LT_24_r$tetrahymena, LT_24_r$loxocephallus) 4.3 Does response diversity predict community stability?
Ready to use ‘comm_data’ to calculate correlations between response diversity measures and stability (and perhaps covariance)!
comm_data %>% glimpse()## Rows: 90
## Columns: 15
## $ study_name <chr> "Leary and Petchey 2009", "Leary and Petchey 2009",…
## $ community_name <chr> "community-1", "community-1", "community-10", "comm…
## $ trait_name <chr> "carrying capacity", "intrinsic growth rate", "carr…
## $ environment_name <chr> "temperature", "temperature", "temperature", "tempe…
## $ cor_respdiv <dbl> 0.5185491, -0.3590885, -0.1281190, 0.9679919, -0.45…
## $ deriv_nss_respdiv <dbl> 1.032674, 1.085592, 1.033750, 1.049787, 1.032133, 1…
## $ deriv_ss_respdiv <dbl> 0.00000000, 0.44585207, 0.09247253, 0.00000000, 0.1…
## $ cv <dbl> 0.4633953, 0.4633953, 0.3565395, 0.3565395, 0.39328…
## $ cov <dbl> 4.104023e-04, 4.104023e-04, 1.277797e-05, 1.277797e…
## $ cor <dbl> 0.71524031, 0.71524031, 0.18853748, 0.18853748, 0.0…
## $ eve <dbl> 0.7642407, 0.7642407, 0.7249706, 0.7249706, 0.59058…
## $ var <dbl> 2.485574e-03, 2.485574e-03, 4.413172e-04, 4.413172e…
## $ total_biomass <dbl> 0.06204328, 0.06204328, 0.03030133, 0.03030133, 0.0…
## $ composition <chr> "colpidium-loxocephallus", "colpidium-loxocephallus…
## $ environment_values <chr> "18-20-22", "18-20-22", "18-20-22", "18-20-22", "18…Simple correlations between CV (and COV) and response diversity (measured as the pairwise correlation, similar to in Leary & Petchey 2009), the dissimilarity (from this study), and the divergence (from this study). I thought the best way to compare would be to recreate the panels with the pairwise responses from Leary & Petchey (e.g., Figs. 4a, 4b, 4f, 4g) and compare with the metrics we propose. I used a simple pearson correlation, as I couldn’t figure out exactly how the correlations were made in Leary & Petchey, but probably some fancier model might be better (maybe even nonlinear fits…). I included the result of the correlation as a footnote in each figure, labelled “Cor = X” where X is the Pearson Correlation coefficient.
# instrinsic growth rate (r) - CV -----
comm_data_r<-comm_data %>%
filter(trait_name %in% "intrinsic growth rate")
panel_A<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate correlation",y = "CV of population biomass",
tag = "a)",caption = str_c("Cor = ",cor(comm_data_r$cor_respdiv,comm_data_r$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = cor_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = cor_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_B<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate (dis)similarity",y = "CV of population biomass",
tag = "b)",caption = str_c("Cor = ",cor(comm_data_r$deriv_nss_respdiv,comm_data_r$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_nss_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_nss_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_C<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate divergence",y = "CV of population biomass",
tag = "c)",caption = str_c("Cor = ",cor(comm_data_r$deriv_ss_respdiv,comm_data_r$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_ss_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_ss_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
# carrying capacity (K) - CV -----
comm_data_K<-comm_data %>%
filter(trait_name %in% "carrying capacity")
panel_D<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity correlation",y = "CV of population biomass",
tag = "d)",caption = str_c("Cor = ",cor(comm_data_K$cor_respdiv,comm_data_K$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = cor_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = cor_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_E<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity (dis)similarity",y = "CV of population biomass",
tag = "e)",caption = str_c("Cor = ",cor(comm_data_K$deriv_nss_respdiv,comm_data_K$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_nss_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_nss_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_F<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity divergence",y = "CV of population biomass",
tag = "f)",caption = str_c("Cor = ",cor(comm_data_K$deriv_ss_respdiv,comm_data_K$cv) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_ss_respdiv, y=cv),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_ss_respdiv, y=cv,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)])
(panel_A + panel_B + panel_C)/(panel_D + panel_E + panel_F)
# instrinsic growth rate (r) - COV ---------
comm_data_r<-comm_data %>%
filter(trait_name %in% "intrinsic growth rate")
panel_A<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate correlation",y = "Summed covariances",
tag = "a)",caption = str_c("Cor = ",cor(comm_data_r$cor_respdiv,comm_data_r$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = cor_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = cor_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_B<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate (dis)similarity",y = "Summed covariances",
tag = "b)",caption = str_c("Cor = ",cor(comm_data_r$deriv_nss_respdiv,comm_data_r$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_nss_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_nss_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_C<-comm_data_r %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Intrinsic growth rate divergence",y = "Summed covariances",
tag = "c)",caption = str_c("Cor = ",cor(comm_data_r$deriv_ss_respdiv,comm_data_r$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_ss_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_ss_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
# carrying capacity (K) - COV -----------
comm_data_K<-comm_data %>%
filter(trait_name %in% "carrying capacity")
panel_D<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity correlation",y = "Summed covariances",
tag = "d)",caption = str_c("Cor = ",cor(comm_data_K$cor_respdiv,comm_data_K$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = cor_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = cor_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_E<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity (dis)similarity",y = "Summed covariances",
tag = "e)",caption = str_c("Cor = ",cor(comm_data_K$deriv_nss_respdiv,comm_data_K$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_nss_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_nss_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)]) +
theme(legend.position = "none")
panel_F<-comm_data_K %>%
ggplot() +
theme_classic(base_size = 14) +
labs(x = "Carrying capacity divergence",y = "Summed covariances",
tag = "f)",caption = str_c("Cor = ",cor(comm_data_K$deriv_ss_respdiv,comm_data_K$cov) %>% round(digits = 3))
) +
geom_smooth(aes(x = deriv_ss_respdiv, y=cov),
method = 'lm',col = 'black',alpha = 0.3) +
geom_point(aes(x = deriv_ss_respdiv, y=cov,
col = environment_values, shape = composition),size = 2, alpha = 0.75) +
scale_shape_manual(values = c(1:5)) +
scale_color_manual(values = Isfahan1[[1]][c(6,4,2)])
(panel_A + panel_B + panel_C)/(panel_D + panel_E + panel_F)
Next we calculate the slope (with confidence intervals) of the relationship between community stability and the various measures of response diversity and display the results in a coefficients plot.
comm_data_long <- comm_data %>%
## make long, so that models of all combinations of
## response and explanatory models can be easily run
pivot_longer(names_to = "resp_div_name", values_to = "resp_div_value",
5:7) %>%
pivot_longer(names_to = "community_variable_name", values_to = "community_variable_value",
5:10) %>%
## scale the explanatory variables so model coefficients are more comparable
group_by(study_name, resp_div_name) %>%
mutate(resp_div_value_scaled = scale(resp_div_value)) %>%
ungroup()
mods <- comm_data_long %>%
nest_by(study_name, trait_name, environment_name,
resp_div_name, community_variable_name) %>%
mutate(models = list(broom::tidy(lm(community_variable_value ~ resp_div_value_scaled,
data = data), conf.int = TRUE))) %>%
select(-data) %>%
unnest(cols="models")
mods$resp_div_name[mods$resp_div_name %in% 'cor_respdiv']<-"Correlation"
mods$resp_div_name[mods$resp_div_name %in% 'deriv_nss_respdiv']<-"(Dis)similarity"
mods$resp_div_name[mods$resp_div_name %in% 'deriv_ss_respdiv']<-"Divergence"
mods$resp_div_name <- mods$resp_div_name %>% as.character() %>% parse_factor(levels = c("Correlation","(Dis)similarity","Divergence"))
mods$trait_name <- mods$trait_name %>% as.character() %>% parse_factor(levels = c("intrinsic growth rate","carrying capacity"))
study_oi <- "Leary and Petchey 2009"
coef_plot <- mods %>%
filter(study_name == study_oi,
term == "resp_div_value_scaled",
community_variable_name == "cv") %>%
ggplot(aes(x = resp_div_name, y = estimate)) +
theme_classic(base_size = 14) +
geom_hline(yintercept = 0,lty = 2) +
geom_point() +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high),fatten = 5)+
labs(title = paste0(study_oi)) +
xlab("Measure of response diversity") +
ylab("Relationship with community stability \n(model slope)") +
facet_grid(cols = vars(trait_name))
coef_plot
# same for COV
coef_plot2 <- mods %>%
filter(study_name == study_oi,
term == "resp_div_value_scaled",
community_variable_name == "cov") %>%
ggplot(aes(x = resp_div_name, y = estimate)) +
theme_classic(base_size = 14) +
geom_hline(yintercept = 0,lty = 2) +
geom_point() +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high),fatten = 5)+
labs(title = paste0(study_oi)) +
xlab("Measure of response diversity") +
ylab("Relationship with community stability \n(model slope)") +
facet_grid(cols = vars(trait_name))
coef_plot2
The model (linear regression) is rather crude (a more sophisticated could account for some of the dependencies among observations). Nevertheless, the point is to show that this approach can be used to assess the performance of different measures of response diversity, and different traits used as the basis of the measure. In this example, response diversity calculated with the intrinsic growth rate trait and as the sign sensitive amount of variation in the derivative has the strongest relationship with community stability, and is the only relationship with confidence interval not overlapping zero.
And now with the correlation
study_oi <- "Leary and Petchey 2009"
mods <- comm_data_long %>%
nest_by(study_name, trait_name, environment_name,
resp_div_name, community_variable_name) %>%
mutate(models = list(broom::glance(lm(community_variable_value ~ resp_div_value_scaled,
data = data)))) %>%
select(-data) %>%
unnest(cols="models")
mods$resp_div_name[mods$resp_div_name %in% 'cor_respdiv']<-"Correlation"
mods$resp_div_name[mods$resp_div_name %in% 'deriv_nss_respdiv']<-"(Dis)similarity"
mods$resp_div_name[mods$resp_div_name %in% 'deriv_ss_respdiv']<-"Divergence"
mods$resp_div_name <- mods$resp_div_name %>% as.character() %>% parse_factor(levels = c("Correlation","(Dis)similarity","Divergence"))
mods$trait_name <- mods$trait_name %>% as.character() %>% parse_factor(levels = c("intrinsic growth rate","carrying capacity"))
coef_plot <- mods %>%
filter(study_name == study_oi,
#term == "resp_div_value_scaled",
community_variable_name == "cv") %>%
ggplot(aes(x = resp_div_name, y = r.squared)) +
theme_classic(base_size = 14) +
geom_hline(yintercept = 0,lty = 2) +
geom_point() +
#geom_pointrange(aes(ymin = conf.low, ymax = conf.high),fatten = 5)+
labs(title = paste0(study_oi)) +
xlab("Measure of response diversity") +
ylab("Relationship with community stability \n(r-squared)") +
facet_grid(cols = vars(trait_name))
coef_plot
# same for COV
coef_plot2 <- mods %>%
filter(study_name == study_oi,
#term == "resp_div_value_scaled",
community_variable_name == "cov") %>%
ggplot(aes(x = resp_div_name, y = r.squared)) +
theme_classic(base_size = 14) +
geom_hline(yintercept = 0,lty = 2) +
geom_point() +
#geom_pointrange(aes(ymin = conf.low, ymax = conf.high),fatten = 5)+
labs(title = paste0(study_oi)) +
xlab("Measure of response diversity") +
ylab("Relationship with community stability \n(r-squared)") +
facet_grid(cols = vars(trait_name))
coef_plot2
5 Results of the systematic review
Here are the results from the systematic review conducted by Samuel RP-J Ross in December 2021-January 2022. Columns are as follows:
Date = Date on which paper was included in review. WoS_search = Search string for reproducing systematic review (NA indicates papers not returned in WoS results). Handled_by = Intitials of author handling each paper during screening (SRP-JR = Samuel RP-J Ross). Included = Binary, included (1) or not (0) as an empirical paper that measures response diversity after all screening. Authors = Short form author list of each article. Title = Title of article. Journal = Journal of article publication. Year = Year of publication. Vol = Volume of publication. Issue = Issue of publication. Pg_start = First page number of publication. Pg_end = Last page number of publication. Article_ID = Identifier of articles without page numbers. DOI = Digital object identifier (DOI) of article. Empirical_study = Does the study use empirical data? 1 = yes, 0 = no, 0.5 = sort of (Theory papers etc.). Measures_RD = Does the empirical (or theoretical) study measure response diversity? 1 = yes, 0 = no, 0.5 = sort of. RD_traits = For papers that measure response diversity, do they measure functional response trait diversity? 0 = no, 0.5 = sort of, FDis = yes, using Functional Dispersion, RaoQ = yes, using Rao’s Quadratic dissimilarity. RD_Int_term = For papers that measure response diversity, do they measure an interaction term between species responding to an environmental variable? 0 = no, 0.5 = sort of, 1 = yes. RD_Other = For papers that measure response diversity, do they take another approach (i.e. not traits or interaction terms)? 0 = no, 1 = yes. Method_Description = Short written description to help clarify method (see individual articles for more details). Notes = Any remaining comments from SRP-JR to help justify inclusion/exclusion, or additional comments.
review <- read_csv(here("data/WoS_results.csv"))
review## # A tibble: 187 × 21
## Date WoS_search Handled_by Included Authors Title Journal Year Vol Issue
## <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <dbl> <dbl> <chr>
## 1 17/12… response … SRP-JR 1 Lalibe… Land… ECOLOG… 2010 13 1
## 2 17/12… response … SRP-JR 1 Craven… Eval… DIVERS… 2016 22 5
## 3 17/12… response … SRP-JR 1 Sasaki… Spec… JOURNA… 2019 107 4
## 4 17/12… response … SRP-JR 1 Morel … Pass… ECOLOG… 2020 30 1
## 5 17/12… response … SRP-JR 1 Aquilu… Eval… ECOLOG… 2020 30 5
## 6 17/12… response … SRP-JR 1 Thornh… The … GLOBAL… 2018 24 7
## 7 17/12… response … SRP-JR 1 Mumme … Func… BIOLOG… 2015 191 <NA>
## 8 17/12… response … SRP-JR 1 Mina e… Netw… ECOLOG… 2021 31 1
## 9 17/12… response … SRP-JR 1 Dobert… Logg… JOURNA… 2017 105 5
## 10 17/12… response … SRP-JR 1 Altoma… Asse… ACTA O… 2021 111 <NA>
## # … with 177 more rows, and 11 more variables: Pg_start <chr>, Pg_end <chr>,
## # Article_ID <chr>, DOI <chr>, Empirical_study <dbl>, Measures_RD <dbl>,
## # RD_traits <chr>, RD_Int_term <dbl>, RD_Other <dbl>,
## # Method_Description <chr>, Notes <chr>